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ABSTRACT 

We describe hydrodynamical simulations of galaxy formation in a A cold dark matter (CDM) cosmology 
performed using a subresolution model for star formation and feedback in a multiphase interstellar medium 
(ISM). In particular, we demonstrate the formation of a well-resolved disk galaxy. The surface brightness 
profile of the galaxy is exponential, with a fi-band central surface brightness of 21.0 mag arcsec" 2 and a scale- 
length of Rd = 2.0/r'kpc. We find no evidence for a significant bulge component. The simulated galaxy falls 
within the /-band Tully-Fisher relation, with an absolute magnitude of / = -21.2 and a peak stellar rotation 
velocity of V rot = 121.3 kms -1 . While the total specific angular momentum of the stars in the galaxy agrees 
with observations, the angular momentum in the inner regions appears to be low by a factor of ~ 2. The 
star formation rate of the galaxy peaks at ~ 7 M Q yr" 1 between redshifts z = 2-4, with the mean stellar age 
decreasing from ~ 10 Gyrs in the outer regions of the disk to ~ 7.5 Gyrs in the center, indicating that the disk 
did not simply form inside-out. The stars exhibit a metallicity gradient from 0.7 Z Q at the edge of the disk to 
1 .3 Zq in the center. Using a suite of idealized galaxy formation simulations with different models for the ISM, 
we show that the effective pressure support provided by star formation and feedback in our multiphase model 
is instrumental in allowing the formation of large, stable disk galaxies. If ISM gas is instead modeled with an 
isothermal equation of state, or if star formation is suppressed entirely, growing gaseous disks quickly violate 
the Toomre stability criterion and undergo catastrophic fragmentation. 
Subject headings: galaxies: evolution - galaxies: formation - methods: numerical 



1. INTRODUCTION 

Numerical simulations have become one of the most impor- 
tant theoretical tools for exploring the complicated problem of 
galaxy formation. M odern numerical work on galaxy forma- 
tion w as initiated b y iKatz & Gunnl (119911) . [Navarro & Benzl 
d!99 lh . and iKatzl i 19921) . who included radiative cooling 
by hydrogen and helium, and attempted to account for star 
formation and feedback processes. Subsequent work em- 
ploy ed_jnore_sophisticated treatments of supernova feed- 
back (Navarro & White 1993), models for chemical enrich- 
ment ( Steinmetz & Miiller 1994), and methods to 'zoom 
in' on galaxies of inter est within cosmological simulations 
JNavarro & White! 19941) 

While these studies achieved limited success in creating ro- 
tationally supported disks, the simulated galaxies typically 
failed to reproduce their observed counterparts. In general, 
the simulated disks were found to be too small an d were more 
centrally concentrated than actu al galaxies ( Nava rro" & White| 
11993111994 iNavarro et al1ll995l) . In addition, the star forma- 
tion in the models was overly efficient, convert ing too la rge a 
fraction of the gas into stars by the pres ent day ( Navarro et al. 
ITo^lSteinrnetz & Mull ed! 994 17995) . 

Many of the shortcomings of the early modeling can be tied 
to the coarse resolution of the simulations and the manner in 
which feedback was treated. In a cosmological context, the 
multiphase structure of the ISM within galaxies cannot be re- 
solved in detail, so star-forming gas is generally described as 
a single-phase medium. If feedback energy is deposited into 
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this gas purely as thermal energy, it will be quickly radiated 
away, and any impact of feedback from star formation will 
be lost. In addition, if star formation is very efficient in low 
mass halos at high redshifts, a large number of dense, compact 
galaxies will form and eventually be incorporated into the in- 
ner regions of larger galaxies through hierarchical merging. 
If these baryonic clumps lose angular momentum to dark ha- 
los, the resulting objects will be smaller in extent than if the 
baryons had been accreted smoothly. Thus, an early collapse 
of baryons opens a channel for gas to lose angular momen- 
tum, yielding galaxies that are too compact, lack angular mo- 
mentum, and produce stars overly efficiently compared with 
observations. 

Subsequent work attempted to alleviate these shortcom- 
ings by employing stronger forms of feedback and modi- 
fications to the cosmological framework. In some cases, 
the problems noted above were exacerbated. For example, 
photoheating by a diffuse ultraviolet (UV) background was 
found to further reduce the angu lar momentum co ntent of 
the simulated galaxies (Navarro & Steinmetzl ll99'7l) . Other 
efforts included the imp act of preheating an d gas blow-out 
from small halos ( Somm er-Larsen et alJ 1999) and the forma- 
tion of galaxies in a warm dark matter (WDM) cosmology 
(Sommer-Larsen & Dolgov 2001). More direct comparisons 
to observations were made possible by incor porating spec- 
tral s ynthesis techniques into the modeling (Contar do et al.l 
1998) and by using the Tully-Fisher relation to constrain the 
hierarchical origin of galax ies dSteinmetz & Navarro! 1 1999t 
iNavarro & Stein metz 2000). However, the galaxies in these 
simulations were still too concentrated. 

The most recent studies of disk formation have yielded 
somewhat more promising results. Using a model of self- 
propagating star formation co mbined with supernov a feed- 
back and a UV background, Sommer- Larsen et alJ {2002, 
120031) produced disk galaxies deficient in angular momentum 
by less than an order of magnitude. iGovernato et alJ (120021) 
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report the formation of realistic disk galaxies in ACDM and 
AWDM cosmologies using simulations that include standard 
prescriptions for cooling, a UV background, and star forma- 
tion. l^ad^^O <EQQl3E) present detailed analyses of simu- 
lated galaxies with kinematic and photometric properties sim- 
ilar to observed Sab galaxies. 

Although the various attempts to simulate disk forma- 
tion have provided some impressive success es, the es sential 
physic s behind this process remains unclear. Governat o et alJ 
(2002) claim that the low angular momentum of sim- 
ulated galaxies in previous works owed at least partly 
to inadequate resolution. In addition, they emphasize 
the impact that the matter power spectrum can have on 
galaxy formation through simulations of WDM universes 
in which the bulge and spheroid components of galax- 
ies are smaller t han in CDM models . These fi ndings 
are supported by [S ommer-Larsen et al. (2002, 2003J) and 
Sommer-Larsen & Dolgov (2001), whose calculations indi- 
cate that high resolution, strong stellar feedback, and warm 
dark matter can p roduce realistic disk galaxies. Finally, 
lAbadi et"aT]|l200 3a) suggest that a crucial ingredient for form- 
ing proper disks is an implementation of feedback that can 
regulate star formation. 

Here, we use a 'multiphase model' to describe star-forming 
gas, in which the ISM consists of col d clouds in pressure 
equilibrium with an ambient hot phase ( Springel & Hernquist 
2003a). Radiative cooling of gas leads to the growth of 
clouds, which in turn host the material that fuels star for- 
mation. The supernovae associated with star formation pro- 
vide feedback by heating the ambient medium and evapo- 
rating cold clouds. The feedback treatment establishes a 
self-regulated cycle for star formation, pressurizing the star- 
forming gas. Implemented in a subresolution manner, our ap- 
proach makes it possible to obtain numerically converged re- 
sults for the star formation rate at moderate resolution. This 
feature is particularly important for simulations of CDM cos- 
mologies, where early generations of galaxies may be difficult 
to resolve. 

Using this model of star-forming gas, Springel & Hernauistl 
(2003b) obtained a converged prediction for the his- 
tory of cosmic star formation that agrees well with 
observations at redshifts z < 4 {Springel & Hernquist 
2003b; Hernquist & Scringe] 2003f). 1 In subsequent work, 
Nagamine et al. (2003a b d) showed that this model also ac- 
counts for the observed abundance and star formation rates of 
damped Lyman-alpha absorbers and Lyman-break galaxies at 
Z~3. 

In what follows, we study the consequences of our multi- 
phase model for the ISM on the formation of disk galaxies. 
First, we employ a high-resolution simulation to identify re- 
alistic disk galaxies in a cosmological context. For one such 
disk galaxy, we examine its kinematic and photometric prop- 
erties in detail, demonstrating good agreement with observa- 
tions of local spirals. Second, we use a set of idealized simula- 
tions to study disk formation in individual dark matter halos to 
isolate the effect of different models for the equation of state 
of the ISM. This analysis demonstrates that the feedback in 
our multiphase model alters the dynamics by pressurizing the 

1 We note an error in figure 12 of Springel & Hernauist 1 2003b) where the 
observational estimates of the star formation rate are plotted too high by a fac- 
tor of h = 1.4. When corrected, the observed points are in better agreement 
with th e theoretical estimates; see astro-ph/0206395 and Nagamine et al. 
120033) . 



star-forming gas and stabilizing forming disks against frag- 
mentation. We emphasize that this aspect of our modeling 
does not depend on the details of our prescription for star 
formation and feedback, but is determined by the effective 
equation of state for star-forming gas. Thus, our conclusions 
should obtain generally, provided that the actual bulk equa- 
tion of state for the ISM has characteristics similar to those of 
our description. 

In § |3 we present our simulation method. We review our 
analysis procedure in § [3] and discuss our findings for the 
structural (§ |4} and kinematic properties (§ [5jl of one simu- 
lated disk galaxy. In § |6] we describe our idealized simu- 
lations and their results. Finally, we conclude and suggest 
directions for further research in §0 

2. COSMOLOGICAL SIMULATION 

Our simulations were performed using the parallel 
iV-body/smoothed particle hydrodynamics (SPH) code 
GADGET2 in its "conser vative entropy" formulation 
( Springel & Hernquist 2002), to mitigate problems with 
lack of ener gy and entropy conservation in older treatments 
of SPH (e.g. lHernquislfTggllO'Shea et alJ l2003). We adopt 
a flat ACDM cosmology with cosmological parameters 
fl m = 0.3, tt A = 0.7, tt h = 0.04, and <j 8 = 0.9, and set the pri- 
mordial power spectrum index to n = 1 . Throughout, we select 
a value for the Hubble constant of Ho = 100/ikms~ 1 Mpc~ 1 
with h = 0.7. 

For our cosmological simulation, we populate a periodic 
volume of 10 /r'Mpc on a side with 144 3 low-resolution 
dark matter (LRDM) particles. At the center of the box, a 
5 /r'Mpc cubic region is selected as a high resolution re- 
gion where we replace the LRDM particles with particles 
of eight-times lower mass. The initial displacement field is 
then calculated following a standard "zooming" procedure 
llTormenl 119971 iPower et all 12003). where small scale per- 
turbations are added appropriately in the high-resolution re- 
gion. Note that the high-resolution zone does not target 
a particular object, with the intent of eliminating bias that 
could be introduced by selecting halos that may be intrinsi- 
cally favorable for disk galaxy formation. We further split 
the high-resolution particles into dark matter (HRDM) and 
gas. The resulting particle masses of each component are 
then m LRDM = 2.79 x 1O 7 /! _1 M , m H R D M = 3.02 x lO^r'M©, 
and m gas = 4.65 x 10 5 /i _1 Mq. We set the gravitational soften- 
ing length for the high-resolution particles to 0.65 comoving 
/r'kpc. While our 10 h~ l Mpc box is too small to be fully 
representative of the z = universe, the volume is sufficient 
for our purposes as our current work does not concern, for ex- 
ample, large-scale correlations of galaxies or the global mass 
function. Here, we are interested only in individual, galactic- 
sized objects which are mostly unaffected by the simulation 
box size. 

We include a UV backgrou nd and radia tive coolin g and 
heatin g in the manner of iKatz et alJ {l996) and [Dave et alJ 
(1999), as well as star formation, supernova feedback, and 
metal enrichme nt. We employ the mu ltiphase model devel- 
oped by Snring eT& Hernqu ist! (|2003al) to describe the star- 
formi ng gas ("see also lYepes et alJl997UHultman & Pharasvn 
1999). Our approach accounts for s ome of the key aspects of 
the multiphase structure of the ISM ( iMcKee & Ostrikerl 19771) 
without spatially resolving the different phases explicitly. In- 
stead, a statistical mixture of the phases is computed analyti- 
cally, taking into account the growth of cold clouds embedded 
in a supernova-heated ambient phase, the formation of stars 
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Table 1 . Springel & Hernquist (2003a) Model Parameters 
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FIG. 1 . — Effective equation of state for star-forming gas in our multiphase 
model i Springel & Hernquist 2003a). The effective pressure, P c ff, is plotted 
in cgs units versus the total gas density in hydrogen atoms per cubic centime- 
ter. The solid line is the exact equation of state, while the thick dashed line is 
a simple fit (see text), which is accurate to about 1%. The vertical dotted line 
shows the transition between an isothermal gas (to the left) and the effective 
equation of state in our multiphase model (to the right). 



out of the cloud material, and the evaporation of clouds in 
supernova remnants. 

Our strategy is motivated by two basic limitations in mod- 
eling star formation on cosmological scales. First, there is no 
fundamental theory of this process. Second, large-scale sim- 
ulations lack the resolution that would be required to char- 
acterize the ISM in detail. Our statistical method decouples 
our ignorance of star formation from its impact on galactic 
scales, by coarse-graining the physics of star-forming gas. In 
principle, this strategy could be applied to any model of star 
formation. 

For a detailed descriptio n of our m ultiphase model we refer 
the reader to Springel & Hernauist (2003a). For definiteness, 
we list in Table 1 the parameter values adopted in our sim- 
ulations. The parameters include the efficiency of the evap- 
oration process Ao, the mass fraction of stars that are short- 
lived and die as supernovae f3, the cold gas cloud temperature 
Tcioud, the effective supernova tempera ture 7sn, and the gas 
consumption time-scale f ( j. As Springel & Hernauist (2003a) 
argue, all of these aside from one, e.g. fj, can be constrained 
by simple physical arguments. We adjust f ( , to match the Ken- 
nicutt Law which describes empi rically observed star forma- 
tion rates in nearby galaxies (e.g. Kennicutt 1989, 1998). 

On the surface, our description of star-forming gas appears 
complicated, with the properties of the different phases de- 
scribed by differential equations that include the various pro- 
cesses listed above. However, in the end, our model can be 
reduced to t wo main ingredients that make our results easier 
to interpret. Sr>rin geT& Hernauist! (I2003al) showed that the 
star-forming gas quickly establishes a self-regulated cycle in 
which the mass and energy of the phases can be approximated 
by equilibrium solutions. In this limit, the impact of star for- 
mation and feedback can be reduced to: 1) the rate at which 
gas is converted into stars, and 2) the effective equation of 
state for star- forming gas. 

Here, as in lSpringel & HernauistH200'3ah . we choose to pa- 



Parameter 


Symbol 


Value 


Supernova Evaporation 


Ao 


1000 


Mass Fraction of Stars > 8Mq 




0.1 


Cold Cloud Temperature 


^clouds 


1000 K 


Effective Supernova Temperature 




10 s K 



rameterize the star formation rate by 

^ = (1-/3A 
df f» 



where 



-1/2 



(1) 



(2) 



The free parameter fg is chosen to match the Kennicutt Law, 
p t h is a threshold density above which gas is subject to ther- 
mal instability and star formation, and p c is the density of 
clouds which, in eq uilibrium, is given by equation (18) of 
Snrin gd"& Hernquist! fc003al) . ISpringel & Hernquist! |2003a) 
showed that p± can be fixed by e.g. requiring the equation of 
state to be continuous. 

In our subresolution model, feedback from supernovae adds 
thermal energy to the ISM, pressurizing it and modifying its 
equation of state. On the scales relevant to the dynamics of the 
gas, this is described by an effective equation of state, which 
can be written 

P e & = (7-l)P"eff, (3) 

where 7 is the ratio of specific heats, p is the total gas den- 
sity, and M e ff is the effective specific thermal energy of the 
star-fo rming gas which, in eq uilibrium, is given by equation 
(19) of Springel & Hernauist (2003a). The effective equation 
of state for the particular model parameters employed in our 
study is illustrated in Figure ^ m c g s units. Also shown in 
Figure^is a simple fit to the equation of state, which is accu- 
rate to about 1 % and is given by 

logoff = 0.050 (logw H ) 3 - 0.246 (logn H ) 2 + (4) 
1.749 logn H - 10.6, for log n H > -0.89 . 

In the particular example shown in Figure ^ the gas is as- 
sumed to be isothermal for densities p < p t ^. In cgs units, 
for the parameter choices adopted here, the critical density at 
which the gas begins to depart from isothermality is «H,th = 
0.128 cm" 3 . This choice ensures that the effective pressure 
is a continuous function of density. For the equation of state 
shown in Figure 1, a different value of nn,th would introduce 
an unphysical jump in pressure at the transition density where 
the gas is no longer isothermal. 

Our implementation of fe edback differs signi fi cantly 
from previous work s (e.g. ISommer-Larsen et ail l2003t 
Governato et al. 2002). In our subresolution model, feedback 
energy is "stored" in the surrounding gas, adding pressure 
support to it. As shown in Figure ^ at densities p > p t h, the 
equation of state becomes stiffer than isothermal due to this 
feedback energy. This pressurization regulates star formation 
occurring within the ISM and, as we discuss later, modifies 
the bulk dynamics of the gas by altering the pressure gradient 
term in Euler's equation. It is this property of our descrip- 
tion that we believe is responsible for many of the differences 
between our numerical results and earlier simulations. 
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FIG. 2. — Entire simulation volume at z = projected onto a 10 h Mpc square area, showing the end-state of the simulation. The upper left and right panels 
show the HRDM and LRDM particles, respectively. The lower left and right panels show the gas and star particles, respectively. Owing to gravitational effects 
the low- and high-resolution particles have intermingled at the boundary of the high-resolution region, but the LRDM particles have not penetrated the inner 
high-resolution region. Star particles have spawned in regions where the gas density has exceeded the critical density for star formation. 



In passing, we note that for the purposes of large-scale dy- 
namics, our multiphase model can be reduced to the choice of 
star formation rate, e.g. equation 0, and the effective equa- 
tion of state of star-forming gas, e.g. equation l|5} and Fig- 
ure [2 The results presented here could be obtained with any 
hydrodynamical code that uses the same expression for the 
star formation rate and an equation of state like that in Fig- 
ure^ without specific reference to all the complexities of our 
multiphase model. 

In Figure |2] we show the evolved simulation volume at 
Z = 0, in projection. The upper and lower left panels show 
the high-resolution dark matter (HRDM) and gas particles, 
respectively. The high-resolution region of the simulation has 
increased in size from the initial 5 h Mpc cubic region at 
the center of the simulation, owing to the growth of gravita- 
tional structure. Some of the low- and high-resolution parti- 
cles have mixed in the boundary region of the low-resolution 
zone, whose shape has also been distorted. However, the 
LRDM particles (upper right panel) have not penetrated into 
the inner regions of the high-resolution region where many 
luminous galaxies have formed out of the gas, as marked by 
the star particles (lower right panel). 

3. DISK GALAXY SELECTION 

We use a standard friends-of-friends (FOF) group-finding 
algorithm to identify dark matter halos, restricting ourselves 



to objects comprised of high-resolution particles. We then 
identify the most bound star or gas particle within each ob- 
ject and produce a new collection of particles for each halo 
that encompasses all particles that lie within the radius /?200 
that encloses 200 times the critical density around these most- 
bound baryonic particles. The group catalogue produced from 
this procedure serves as our primary galaxy sample. 

From our galaxy catalogue, we have visually inspected the 
10 largest objects, measured their angular momenta, and se- 
lected a large disk galaxy for further detailed analysis in this 
paper. This galaxy displays the most disk-like morphology of 
all the galaxies in the simulation resolved with a similar num- 
ber of particles. While we do find other flattened systems in 
the simulation, the majority of them do not display such well 
defined disk-like morphology and kinematics as the galaxy 
described below. 

The galaxy analyzed here, which we term "galaxy CI" for 
simplicity, has 61,372 dark matter, 21,506 gas, and 88,138 
star particles within /?2oo at z = 0. Note that the star particle 
mass was set to half the original mass of the gas particles, 
so the overall baryon fraction within R200 is 14.1%; slightly 
larger than the universal baryon fraction of 13.3%. Figure 
[5] shows the star (left panels) and gas (right panels) particle 
distributions within 0.1 R200 at z = 0, projected perpendicular 
(upper panels) and parallel (lower panels) to the direction of 
the stellar angular momentum. In addition, the lower right 
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FIG. 3. — Disk galaxy at z = 0, displaying the star (left panels) and gas (right panels) particle distributions within O.IR200. projected perpendicular (upper panels) 
and parallel (lower panels) to the direction of the stellar angular momentum vector. The lower right panel also traces the gas velocity field and demonstrates the 
rotational support of the galaxy. 



panel traces the gas velocity field in the plane of the disk. The 
gas in the galaxy has collapsed into a rotationally supported 
disk, surrounded by a thicker stellar disk. The velocity field of 
the galaxy illustrates the mainly circular trajectories of the gas 
particles in the disk plane. A weak bar is visible in the disk 
gas, but is not strong in the stars that dominate the baryonic 
mass of the galaxy. We note that the stellar component is 
slightly elongated in the direction of the bar, however. The 
galaxy has reached a mass of M 2 oo = 2.16 x 10 n h~ l M Q at 
z = 0, with a stellar mass of M* = 2.05 x I0 w h~ l M© and a 
gas mass of M gas = 9.44 x 10 9 A -1 M©. The baryonic mass of 
the galaxy is then 68% stellar and 32% gaseous. 

4. STRUCTURAL PROPERTIES OF GALAXY CI 

4.1. Mass surface density profiles 

Figure |4] shows the stellar (diamonds) and gas (triangles) 
mass surface densities of galaxy CI projected along the stellar 
angular momentum vector onto the disk plane. The stellar 
profile includes particles within 30/i _I kpc of the disk center, 
while the gas profile includes particles only within 2 h~ l kpc 
of the disk plane. The stellar and gas mass surface densities 
are roughly exponential out to r = 8/i _I kpc. If we fit them 
with an exponential of the form 

£(r) = £oe- r / R ", (5) 
we find best fit values of Eo,* = 452 h M Q pc" 2 , * = 
2.3/r'kpc for the stars and So. gas = 15.3/z M©pc~ 2 , Bd, gas = 



4.7/r'kpc for the gas. In the plot, we include least squares 
fits to the stellar and gas mass surface densities. For com- 
parison, the stellar disk of the Milky Way has a scale-length 
Smw ^3.5 kpc, similar to the physical scale-length of ~ 3.3 
kpc for our simulated galaxy CI. We note that the gas mass 
surface density begins to diverge from the exponential fit be- 
tween 8 — 10/T 1 kpc, at the edge of the gaseous disk. 

4.2. Photometric properties 

From the known masses, metallicities, and ages of the star 
particles, we use a population spectral synthesis code to calcu- 
late the SDSS ugriz JFukugita et alJl996T) AB system magni- 
tudes and JHK magnitud es for each par ticle. The population 
synthesis code assumes a Kroupa ( 2001 ) IMF and limits input 
metallicities to a range of 0.005 Z© < Z < 2.5 Z©, outside of 
which we adop t the m inimum or maximum value. We then 
use the Fukugit aet al.l (fl996) conversion s from SDSS colors 
to the Johnson-Morgan-Cousins UBVRI dJohnson & M organ 
1953; Cousins 1978) Vega system to calculate the absolute 
magnitude of the star particles in each band. We bin star par- 
ticles within 30 h~ l kpc of the center of the galaxy into an- 
nuli and calculate their total luminosities to determine a sur- 
face brightness profile. The measured surface brightness of 
galaxy CI in B, V, I, and /T-bands is plotted in Figure|5] The 
best fit B-band central surface brightness is /io b = 21.0 mag 
arcsec" 2 with a B-band luminosity scale-length of R c i.b = 2.0 
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FIG. 4. — Mass surface density for the stars (diamonds) within 30/i~' kpc of the galaxy center and gas particles (triangles) within 2/T 1 kpc of the disk plane. 
Using an exponential profile S(r) = So exp(-r/K,;) to fit each component separately, we find best fit values of So,* = 452ft Mqpc -2 , Rd.+ = 2.3/r'kpc for the 
stars and So. 6as = 15.2ft Mqpc~ 2 , S gas = 4.7/r'kpc for the gas. The least squares fit for the stars (solid line) and gas (dashed line) are plotted over the measured 
values. 

FIG. 5. — Surface brightness profile of the galaxy in B (stars), V (dashed line), / (dashed-dotted line), and K (dotted line) bands. The solid line is the fit to the 
B-band profile, with best fit central surface brightness /^o,b = 21.0 mag arcsec -2 and B-band disk scale length Rj B = 2.0ft _1 kpc. The standard central B-band 



surface brightness value for disk galaxies is ^b.Fic 
the galaxy. 



= 21.65 mag arcsec . The disk brightness is exponential in each band throughout the inner 10 h~ l kpc of 



h 1 kpc. Our measured value for the B-band central surface 
brightness agrees well with the canonical value of /xo,b = 2 1 .65 
mag arcsec" 2 for spiral galaxies (Freem anll970l) . 

Figure|5]demonstrates that the surface brightness profile of 
galaxy CI is exponential from the very center of the galaxy 
out to r ~ 10/i -1 kpc, and shows no evidence for a bulge com- 
ponent. To our knowledge, this object is the first published 
example of a galaxy formed within a cosmological simulation 
that displays an exponential surface brightness profile with no 
significant bulge. 

For reference, the u, g, v, r, i, z, U, B, R, I, J, H, and 
^-band absolute magnitudes of the galaxy, measured by sum- 
ming the luminosities of all stars used to calculate the surface 
brightness profiles out to 25.0 mag arcsec -2 in each band, are 
listed in Table 2. The galaxy has realistic colors, and its prop- 
erties agree with recent 2MASS determinations of B-K vs. 
Mk color measurements of nearby spiral galaxies (J. Huchra 
2003, private communication). We also list the stellar mass- 
to-light ratios T x in each band, for comparison with observed 
galaxies. The stellar mass-to-light ratios in the g, r, i, and z 
bands agree well with the properties of observed galaxies in 
the Sloan survey JKauffmann et alJ2003l) . 



formed early on, with the majority of stars in the galaxy at 
z = being older than 9 billion years. A large fraction of these 
stars formed in progenitor halos that later merged to form the 
galaxy. 

Figure illustrates the history of the metal content of the 
galaxy. Each panel shows the average metallicity of newly 
formed stars (solid line) over the lifetime of the galaxy and 
the average metallicity of all stars formed before each epoch 
(dashed line). While the metallicity of newly formed stars 
varies in a complex manner, reflecting the details of gas ac- 
cretion and merging, the average metallicity of the galaxy 
evolves smoothly to the current epoch. The large feature at 
z = 0.8 corresponds to a major merger in the assembly his- 
tory of the galaxy that significantly altered the metallicity be- 
tween z = 0.7 and z = 0.5. Note that stars formed before 9 
billion years ago have average metallicities of < 0.8 Z Q . We 
hence expect to see an older population of relatively metal- 
poor stars, while regions of the galaxy that actively formed 
stars even at late times should feature a younger, metal-rich 
population. If a disk does not form stars at the same rate ev- 
erywhere, we would expect to see stellar age and metallicity 
gradients between regions of differing star formation activity. 



4.3. Star and metal formation histories 

Within our simulation code, we record the formation time 
of each star particle, allowing us to infer the detailed star for- 
mation history of each object. Similarly, the metallicity of the 
star particles can be used to determine the metallicity evolu- 
tion of the galaxies. Figure[6]shows the star formation history 
of galaxy CI. As a function of redshift (left panel), the star 
formation rate rises below redshift z = 10 and peaks at a value 
~ 7 M Q yr" 1 at redshift z ~ 3, after which it declines sharply 
to the present time. In the right panel of Figure[6] we also plot 
the star formation rate against lookback time. Viewing the star 
formation history in this way illustrates that most of the stars 



4.4. Stellar age and metallicity gradients 

Combining the formation time and metallicity of each star 
particle with information about its position allows us to deter- 
mine stellar age and metallicity gradients in the galaxy. The 
left panel of Figure [8] shows the stellar age gradient in the 
disk, measured for all star particles within 15 /T 1 kpc from 
the center of the galaxy. The average age of the star particles 
decreases from ~ 10 Gyr in the outer regions of the disk to 
~ 7.5 Gyr in the center. Beyond a radius of ~ 10/T 1 kpc, the 
average age of the stars is dominated by contributions from 
the stellar halo, where stars are significantly older than in the 
disk. 
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Table 2. Photometric Properties of Galaxy CI 
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FIG. 6. — Star formation history of the entire simulated galaxy as a function of redshift (left panel). The star formation rate in the galaxy increases from 
redshift z = 10 to a peak of ~ 7 Mq yr~' at z ~ 3. The right panel shows the galactic star formation history as a function of lookback time. Most of the stars in 
the galaxy formed earlier than 9 Gyr ago. 
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FIG. 7. — Metal formation history of the simulated galaxy. In each panel, the solid line shows the average mass-weighted metallicity of newly formed stars 
over the lifetime of the galaxy. The dashed line shows the running mass-weighted average metallicity of all stars within the galaxy. The left panel is plotted as a 
function of redshift and the right panel is plotted against lookback time. The large feature at z = 0.8 corresponds to a major merger in the assembly history of the 
galaxy that significantly altered the metallicity between z = 0.7 and z = 0.5. 
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FIG. 8. — Stellar age gradient in the disk (left panel), measured by averaging the ages of star particles in 0.5/i~'kpc radial bins. The average age of the stars 
decreases from ~ 10 Gyr in the outer regions of the disk to ~ 7.5 Gyr in the center. The right panel shows the stellar (solid) and gas (dashed) metallicity gradients 
in the disk. The metallicity of the stars increases from ~ 0.75 Zq at the edge of the disk to ~ 1.3 Zq in the center, while that of the gas in the disk increases 
strongly from primordial values in the outer disk to ~ 2.5 Zq in the central regions. Stars forming after z = 0.25 (dotted line) track the metallicity gradient of 
star-forming gas in the disk. 



In the right panel of Figure [8] we show the metallicity gra- 
dient in the galaxy, plotting the average metallicity of the star 
(solid line) and gas (dashed line) components in 0.5/T 1 kpc 
radial bins. The metallicity gradient of the stars tracks the 
stellar age gradient closely, remaining at a roughly constant 
value of ~ 0.7 Z Q from r = 15h~ l kpc to r = 8/T 1 kpc and 
then increasing to a central value of Z= 1.3 Z Q . Again, the 
halo stars beyond lO/r'kpc are older and deficient in metals 
compared with the younger stars in the disk. Interestingly, 
the gas metallicity gradient is much stronger than the stellar 
metallicity gradient, increasing from close to primordial val- 
ues at r = 8/T 1 kpc to ~ 2.6 Z© near the center of the disk. 
The gas in the inner regions has been enriched by stars form- 
ing at late times (dotted line, for stars forming after z = 0.25), 
whose metallicities have been enhanced by the long history 
of star formation in the galaxy. The metallicity gradients of 
the star-forming gas and recently formed stars in the disk then 
track each other, again highlighting the strong metallicity gra- 
dient in the gaseous disk. 

The active star formation in the disk produces a collection 
of stars that are younger and more metal enriched than in 
the outer regions of the galaxy. Hence, the disk is not sim- 
ply built up inside-out, but i nstead partially forms outside-in 
( Sommer-Larse net alJl2003l) . The more active star formation 
at late times in the inner disk makes this region actually quite 
young, on average. In the halo of the galaxy, where star for- 
mation is dormant, stars are significantly older and more metal 
poor than in the disk. 

5. KINEMATIC PROPERTIES OF GALAXY Cl 

5.1. Angular momentum 

Perhaps the most important kinematic property of disk 
galaxies is their angular momentum, as this quant ity deter- 
mines the sizes of thin, rotationally supported disks (Mo et al. 
1998). While measuring the specific angular momentum of 
each star particle is straightforward, deciding which particles 
to use for measuring the angular momentum content of the 
disk is less clear, in particular when comparing with observa- 



tions. 

We wish to compare our simulat ed galaxy with the /-band 
Tul ly-Fisher relation measured by iMathewson et alJ (119921) 
and Courteau (1997) and the c orresponding j*-y io t relation 
(as compiled by Navarro 1998). The measured rotation veloc- 
ity Vrot used for these relations is determined observationally 
at 2.2 times the exponential scale-length of the stellar surface 
brightness, Rj. We then take V rot to be 



GM(< R v ) 



V m , . 

Ky 

Rv = 2.2Rd.B, 



(6) 
(7) 



where Rj B is the B-band surface brightness scale-length de- 
termined in § 14.21 In what follows, when we refer to V lot we 
always mean the expected rotation velocity corresponding to 
the circular velocity produced by all mass within 2.2 x Rj.b = 
0.045^200 for galaxy Cl. 

Having determined the rotation velocity V rot , we then con- 
sider 3 different approaches for measuring the specific stellar 
angular momentum j+, providing slightly different interpre- 
tations of the angular momentum content of the galaxy. For 
all the star particles, we define the total specific stellar angular 
momentum yVtot to be a sum over the star particles of the form 



7*,tot — 



N, 

^2 m i V i X y i 

1=1 



(8) 



where TV* is the total number of star particles within Rum, v, is 
the velocity with respect to the center of mass of the galaxy, 
r, is the distance from the center of the galaxy, and is the 
mass of the / th star particle. 

We also consider the specific angular momentum ^.r v of 
the stars within the radius Ry at which the rotation velocity is 
measured. This is determined by simply restricting the sums 
in equation (|8j to the stars within a distance Ry of the center 
of the galaxy. 

The final method we employ mimics the observational esti- 
mate for the specific angular momentum content of an expo- 
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nential disk, viz. 



A.obs - 2Rd,B Vic 



(9) 



Unlike the previous two methods, this technique indirectly es- 
timates the angular momentum content by combining a mea- 
surement of the size of the disk with a typical circular veloc- 
ity. While allowing an inter-co mparison of observed galaxies 
with similar rotation velocities ( Abadi et al. 2003a), this ap- 
proach implicitly assume s that the disk has the rotation curve 
of an exponential disk (ICourteaulll997l) . which is not true in 
general. 

We summarize the angular momentum of our simulated 
galaxy in Table 3, listing the values for V rot and the three 
different measurements of j* along with their values in 
physical units. With a rotation velocity V rot = 121.3kms _I 
and total specific stellar angular momentum j+ = 204.2 /T 1 
kpc km s , this galaxy agrees favorably with the observed 
rang e of logf/*) ~ 2.4-3. 2 for galaxies with logfy mf ) ~ 
2.1 JMathewson et afl ll992t ICourteaulH997t iNavarrd 1998). 
However, the region of the galaxy within Ry has a specific 
angular momentum lower by a factor of ~ 2 with j^g v = 
101.0/i~ 1 kpckms~ I . The lower angular momentum of this 
part of the galaxy results from the thickness of the stellar disk. 
Finally, the agreement between the value 7* j0 bs =485 .2/T 1 kpc 
km s" 1 and the observed range of demonstrates that the 
stellar disk has the proper size when compared with observed 
galaxies of similar rotation velocity. 

Overall, we find that the total angular momentum content of 
the galaxy is appropriate for its rotation velocity, when com- 
pared to observations, but that the thick disk reduces the spe- 
cific stellar angular momentum content in the inner regions. 
We stress that the stellar disk of this galaxy should be consid- 
ered deficient in angular momentum content and caution that 
only measuring the total specific stellar angular momentum 
content of a simulated galaxy can lead to an overestimate of 
the specific stellar angular momentum in the central regions 
of the galaxy. 

At present, we cannot determine whether the low angular 
momentum in the inner regions of the galaxy is purely a phys- 
ical effect, or is at least partly a consequence of our relatively 
poor spatial resolution. The gravitational softening length in 
the high-resolution region of our simulation is 0.65/i _I kpc. 
For h = 0.7, this means we cannot resolve thin stellar disks of 
galaxies like the Milky Way. It is possible that our estimate 
of the angular momentum of the inner parts of galaxy CI is 
underestimated because of this problem. 

5.2. Rotation curve 

In Figure [9] we plot the rotation velocity of gas parti- 
cles from within 2 h~ l kpc of the disk plane (crosses), pro- 
jected onto 0.5 h~ l kpc wide annuli. We find that the rota- 
tion velocity of the gas closely traces the total rotation ve- 
locity V(R) = [GM(< R)/R] l l 2 (solid line), as expected for a 
thin disk without a dominant bulge component. If the cen- 



Table 3 . Kinematic Properties of Gal ax y C 1 
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FIG. 9. — Rotation velocity of the gas particles within 2 h 1 kpc of the 
disk plane (crosses), projected onto 0.5 h kpc annuli. We find that the 
gas rotation velocity closely traces the total rotation velocity V(R) = [GM(< 
R)/R] l l 2 (solid line), as expected for a thin disk without a dominant bulge 
component. 

tral regions of the galaxy contained a kinematically important 
bulge, the gas rotation curve would either rise steeply or have 
a central maximum and decline outwards. We interpret the 
slowly rising rotation velocity of the gas as further evidence 
that galaxy CI consists of a rotationally supported exponen- 
tial disk without a significant bulge component. 

5.3. l-band Tully-Fisher relation 

From the photometric modeli ng d escribed in § 14.21 and 
the kinematic measurements in § 15.11 it is possible to com- 
pare galax y CI with the /-band T u lly-Fisher relation as deter- 
mined bvlGiovanelli et all ( 119971) . lHan & Mouldl d!992h . and 
Mathewson et al. II1992I) . For our measured /-band magni- 
tude of / = -21.2 and rotation velocity V 10t = 121.3 km s -1 , 
we find our galaxy to be in good agreement with the aver- 
age observed /-band magnitude of / ~ -20.8 for galaxies with 
logCVk,) ~ 2.1. 

6. ISOLATED HALO SIMULATIONS 

The formation of a well-defined disk galaxy without a dom- 
inant bulge can be viewed as a significant success towards 
obtaining realistic galaxies in cosmological simulations. It is 
therefore of particular interest to understand which aspects of 
our modeling are mainly responsible for this outcome. 

A novel feature of our approach lies in our use of a sub- 
resolution model for the multiphase structure of the ISM, to- 
gether with an effective equation of state that this model im- 
plies for the dense star-forming gas. As can be seen in Figure 
n when the gas is sufficiently dense to be thermally unsta- 
ble and subject to star formation, the equation of state departs 
from that of an isothermal gas, and becomes stiffer with in- 
creasing density. This will affect the bulk hydrodynamics of 
the gas, through the pressure gradient term in the Euler equa- 
tion. Consequently, dense gaseous disks will experience addi- 
tional pressure support from supernova feedback, in a manner 
that is different from previous attempts to model galaxy for- 
mation. 

In order to examine the importance of this effect further, 
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FIG. 10. — Stellar mass distribution, adaptively smoothed and seen edge- and face-on, for the "iso-ism" model at (clockwise from upper left) T = 0.5 Gyr, 
T = 1.0 Gyr, T = 1.5 Gyr, and T = 2.0 Gyr. As the disk develops, gas quickly cools to roug hly 10 4 K and forms stars. The 10 4 K gas is Toomre unstable and 
fragments into dense clumps, which in turn efficiently form stars, leading to the clumpy stellar mass distribution in the disk. As the galaxy ages, the clumps lose 
angular momentum through dynamical friction, coalesce, and fall into the center of the galaxy. 



we use a set of simulations to consider disks forming in iso- 
lated halos (e.g. Spri ngel & Hernauisi20 03a). where we vary 
the description of star-forming gas and associated feedback 
effects. These results show that the physics of dense, star- 
forming gas has a substantial impact on the dynamics of disks, 
indicating that the model of the ISM plays a crucial role in our 
successful formation of a disk galaxy in a cosmological con- 
text. 

6.1. Initial conditions and simulation set 

In order to isolate the dependence of disk formation on 
the physics of the ISM, we investigate the dynamics of gas 
in static spherically symmetric NFW ( Navarro et al J [19951) 
dark matter halos of mass I0 l2 h~ l Mq and concentration 
cnfw = 20. We place 40,000 gas particles that initially trace 
the dark matter in these halos and set the initial temperature 
profile such that the gas is in hydrostatic equilibrium when 
evolved with non-radiative hydrodynamics. We add angular 
momentum to the gas, corresponding to a relatively high spin 

parameter of A = J\E\ 1 ^ 2 /(GM^ 2 ) = 0.1. We then evolve the 
gas forward in time including radiative cooling. Under these 
conditions, we expect the gas to collapse smoothly into a thin, 
rotationally supported disk. We do not intend for this to be 
a realistic model of disk formation in a hierarchical universe, 
but the idealized nature of the simulations described below 
makes it possible to examine the consequences of varying the 



gas physics in detail. 

This set-up provides initial conditions for a suite of simu- 
lations, where we consider three different descriptions of the 
ISM. The first model, which we will refer to as "justcool," 
consists of simply allowing the gas to cool with primordial 
abundances without allowing for any star formation or feed- 
back. The "justcool" scenario serves as an extreme reference 
case where star formation and feedback cannot affect the gas 
dynamics. Since the dense gas that settles into a disk can cool 
efficiently, its temperature will remain close to 10 4 K, effec- 
tively obeying an isothermal equation of state. 

Our second description of the ISM, which we refer to as 
"iso-ism," includes gas cooling and star formation, according 
to equation 0, but no feedback effects. Gas in the disk again 
obeys a nearly isothermal equation of state in this case, but 
compared to our "justcool" model, dense gas is converted into 
collisionless star particles. For ease of comparison, we use the 
same star formation rate for gas of a given density as in our 
multiphase model. 

Finally, our third descr iption of the ISM is the 
Springel & Hernquist (2003a) multiphase ISM model, which 
includes radiative cooling, star formation, and supernova 
feedback in the form of thermal energy input and cloud evap- 
oration, as described earlier. We refer to this model as the 
"multiphase-ism" model. The primary difference between the 
"multiphase-ism" and "iso-ism" scenarios is that in the mul- 
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FIG. 1 1 . — Stellar mass distribution, adaptively smoothed and seen edge- and face-on, for the "multiphase-ism" model at (clockwise from upper left) T = 0.5 
Gyr, T - 1.0 Gyr, T = 1.5 Gyr, and T = 2.0 Gyr. In the "multiphase-ism" run, the gas cools into a thin disk which begins forming stars at a rate regulated by 
supernova feedback. While the gas is comprised of clouds cool enough to form stars, the contributions from the hot surrounding medium to the effective equation 
of state of the gas provides pressure support to the ISM to prevent Toomre instability. The end result is a large disk galaxy with only a modest bulge component. 



tiphase model feedback from supernovae modifies the effec- 
tive equation of state of star-forming gas, providing increased 
pressure support, as indicated by Figure [2 The stiffer equa- 
tion of state helps to regulate star formation so that it occurs 
at a rate consistent with observations of isolated disk galax- 
ies (lKennicutllll989t 119981 iMartin & Kennicuttl kOOl) once 
the free parameter tn of the m odel is properly adjusted (see 
Snri ngel & Hernauisi2003albl) . 

For each ISM model, we ran a number of simulations, vary- 
ing the numerical integration parameters to test the robustness 
of the results. In particular, we varied the number of SPH 
neighbors and the value of the time-step taken by the code. 
None of these choices had a significant effect on the outcome. 
Therefore, our results appear to have converged with respect 
to time integration. When the number of SPH neighbors is 
varied, instabilities in the disks tend to occur at slightly dif- 
ferent times and locations, but the qualitative behavior of the 
simulations remains unchanged. 

6.2. Comparison of disks formed in isolated halos 

The time evolution of simulations with our three ISM mod- 
els reveals interesting differences. In Figure[K)| we show the 
stellar mass distribution, adaptively smoothed and seen edge- 
and face-on, for the "iso-ism" model at (clockwise from upper 
left) T = 0.5 Gyr, T =1.0 Gyr, T =1 .5 Gyr, and T = 2.0 Gyr. 
Gas in the inner parts of the halo quickly cools to roughly 10 4 



K, settles into a thin, rotationally supported disk, and forms 
stars. The 10 4 K gas in the disk is highly compressible, obey- 
ing a nearly isothermal equation of state. Consequently, the 
self-gravitating gaseous disk becomes gravitationally unsta- 
ble as soon as it reaches a moderate surface density. The 
Toomre ( 1964) instability quickly causes a fragmentation of 
the smooth gas disk into clumps, which collapse further under 
their own self-gravity and form stars on a short time-scale. As 
the galaxy ages, the stellar clumps coalesce through dynam- 
ical friction into ever larger lumps, destroying the disk. By 
T = 2 Gyr, the stars are mainly confined in two large blobs, 
with numerous smaller clumps orbiting the center. 

The evolution of the galaxy in the "iso-ism" model contrasts 
sharply with that in the corresponding "multiphase-ism" run. 
We illustrate this in in Figure ^] Again, the gas cools into a 
thin, star-forming disk. However, supernova feedback regu- 
lates the rate of star formation and the dynamics of the dense 
star-forming gas. Owing to feedback, the ISM is pressur- 
ized, stabilizing the gas against the Toomre instability. Con- 
sequently, a smooth distribution of gas is maintained in the 
plane of the disk and stars form steadily, increasing the size 
of the disk from 6/i~'kpc at T = 0.5 Gyr to > 12/i _1 kpc at 
T = 2.0 Gyr. Since the gas remains stable and does not frag- 
ment into clumps, the disk largely avoids angular momentum 
loss and the simulation yields a large disk galaxy with only a 
modest bulge. 
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FIG. 12. — Compaiison of the gas mass distribution in the "justcool" (left), "iso-ism" (center), and "multiphase-ism" (right) simulations at T = 1.0 Gyr. Both 
the "justcool" and "iso-ism" models produce gas disks fragmented by catastrophic Toomre instability. By contrast, the "multiphase-ism" model yields a stable 
disk supported by additional pressure supplied by the Snrineel & Hernauist 1 2003a) model for the ISM. 



In Figure ^] we compare the gas mass distributions in 
the "justcool" (left), "iso-ism" (center), and "multiphase-ism" 
(right) simulations at T = 1.0 Gyr. Both the "justcool" and 
"iso-ism" models produce clumpy gas disks that fragmented 
by the Toomre instability. In fact, the time evolution of the 
"justcool" and "iso-ism" models is very similar, except that 
in the former case the instabilities occur earlier, showing that 
star formation can have a stabilizing effect by reducing the gas 
surface density. Note, however, that only the "multiphase- 
ism" model produces a gaseous disk that exhibits long-term 
stability. 

In summary, our simulations of the formation of disks in 
isolated halos clearly show that the modeling of the ISM has 
a strong effect on disk stability. In particular, an isothermal 
equation of state for star-forming gas produces gaseous disks 
that are violently unstable. Under these conditions it would be 
difficult if not impossible to form large, extended disk galax- 
ies unless the gas surface density is always kept extremely 
low, for example by adopting an unrealistically short star for- 
mation time-scale that would conflict with the observed Ken- 
nicutt law. 

In contrast, the multiphase model for the ISM yields disk 
galaxies that are more stable. The stiffer equation of state pro- 
duces disks that are more resilient against Toomre instability 
than in models with an effectively isothermal ISM. Fragmen- 
tation is avoided, enabling the steady formation of large stellar 
disks even when the disk is gas-rich and has a high gaseous 
surface density. Therefore, we suggest that the behavior of 
the equation of state in our multiphase description of the ISM 
is the primary reason why our cosmological simulation can 
produce disk galaxies with negligible bulge components. This 
conclusion should be insensitive to the detailed features of our 
multiphase model of the ISM, provided that the actual equa- 
tion of state for star-forming gas be haves similarly to the one 
proposed by Spring eT& HernauislH2003al) . 

7. CONCLUSIONS 

We have performed a high-resolution cosmological simula- 
tion that included gas dynamics, radiative heating and cool- 
ing, and star formation and supernova feedback, modeled in 



the framework of a subresolution description of a multiphase 
ISM. Inspecting the ten largest halos, we identified a well- 
resolved disk galaxy at redshift z = 0. The disk galaxy has 
realistic photometric properties, including a reasonable cen- 
tral surface brightness and exhibits good agreement with the 
/-band Tully-Fisher relation. Most important, this object is 
the first published example of a disk galaxy from a simulation 
that features a simple exponential surface brightness profile 
without having a large bulge component. Our work demon- 
strates that extended disk galaxies can form in a ACDM cos- 
mology within full hydrodynamic simulations, an important 
step forward in understanding the origin of disk galaxies. 

While the total angular momentum of the stellar compo- 
nent of the disk galaxy is appropriate for observed galaxies 
with similar rotation velocities, the inner regions of the galaxy 
have low specific stellar angular momentum. It is not clear 
whether this is a problem with the physics, e.g. excessive an- 
gular momentum loss by the gas, or is at least partly an artifact 
of our inability to completely resolve the thin stellar disk. The 
disk galaxy in our simulation does not have an obvious prob- 
lem with its radial size, but statistical studies with samples of 
simulated disk galaxies will be needed to settle this question. 
However, we find the success achieved here encouraging and 
suggest that the problem with disk size is alleviated by em- 
ploying equations of state for star-forming gas that are stiffer 
than isothermal. 

Using a suite of isolated galaxy formation simulations, we 
have explicitly demonstrated the consequences our treatment 
for the ISM has on the formation of disk galaxies. Simulations 
with an effectively isothermal ISM yield gravitationally un- 
stable disks and fail to produce large, smooth stellar disks. In 
contrast, the multiphase model supplies enough pressure sup- 
port to the star-forming gas to prevent catastrophic Toomre 
instability, allowing stable galactic disks to form. We have 
also verified that these conclusions do not depend on details 
of our numerical integration. 

Our work differs in a number of respects from earlier nu- 
merical attempts to study disk galaxy formation. We have 
used a novel formulation of smoothed particle hydrodynamics 
that, by construction, conserves both energy and entropy si- 
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multaneously, even when smoothing lengths vary in response 
to changes in density (Springel & Hernauist 2002). While we 
consider it unlikely that previous numerical work on galaxy 
formation with SPH was strongly affected by inaccuracies 
in e.g. entropy conservation, it is clearly preferable to use a 
fully conservative scheme that mitigates against any numeri- 
cal "overcooling." 

More important, we have formulated the consequences of 
feedback using a subresolution model for star formation in a 
multiphase ISM. In this scheme, the local supernova feedback 
modifies the effective equation of state for star-forming gas, 
leading to an efficient self-regulation of star formation and a 
stabilizing effect on highly overdense gas. This aspect of our 
modeling feeds directly into the dynamics of forming disks 
by modifying the pressure gradient term in the hydrodynamic 
equations of motion. 

Our formulation of the ISM physics has a number of advan- 
tages over earlier treatments of star formation and feedback. 
Other workers have implemented techniques to allow feed- 
back to have a significant dynamical impact by, e.g. deposit- 
ing thermal energy into gas surrounding star-forming regions 
and artificially delaying the radiative loss of this energy by 
the gas for an ad hoc time interval. While similar in spirit 
to our method, these approaches make the comparison of re- 
sults from different calculations challenging because the con- 
sequences of feedback are difficult to summarize in a concise, 
quantitative form. When expressed in terms of an effective 
equation of state, as done here, the dynamical implications of 
feedback become easier to interpret. A discussion in terms of 
the effective equation of state is also relatively insensitive to 
the detailed physics responsible for pressurizing the gas, and 
enables one to relate the description of feedback to observa- 
tions of the ISM. 

In order to demonstrate the sensitivity of galaxy formation 
to the physics of star-formin g gas, we employed a simple 
multiphase model of the ISM (Springel & Hernauist 2003a), 
where the gas is imagined to be thermally unstable at suf- 
ficiently high densities and exists in two distinct phases, in 
pressure equilibrium. In detail, it is not plausible that this is 
a true picture of the ISM of disk galaxies, because this model 
ignores, e.g. turbulent motion, magnetic fields, and cosmic 
rays. However, we believe that the basic points of our mod- 
eling do not depend on these complications. The evolution 
of the different models shown in §[6] depends mainly on the 
effective equation of state for the star-forming gas on scales 
which characterize the Toomre instability. According to the 
Toomre criterion, a thin sheet of gas in differential rotation 
will be gravitationally unstable if 



where c s is the sound speed, k is the epicyclic frequency, and 
£ is the surface density. In terms of an effective temperature 
of the gas, 

q « i.ir 5 1/2 K35 ^, (ii) 

where K35 = K/(35kmsec~ 1 kpc~ 1 ) and £75 = £/(75MqPC~ 2 ) 
are values characteristic of the Milky Way in the solar neigh- 
borhood, and 7g = r/(10 5 K). If the gas obeyed an isothermal 
equation of state with 7g =0.1, then Q ~ 0.3, and the disk 
would be violently unstable, as for example in the case of our 
"iso-ism" simulations described in § 6. With our equation of 
state, however, the effective temperature at densities charac- 
teristic of a forming galaxy would be T5 ~ 1 (see Fig. 1 of 



Springel & Hernauist 2003a), so Q > 1 and the disk would be 
stable, explaining why our "multiphase-ism" models in § 6 
produce realistic disks. Moreover, for thin disks of gas, the 
fastest growing unstable mode has a wavelength which is a 
signifi cant fraction of the disk size (e.g. lBinnev & Tremaind 
1987) and, therefore, the growth of the instability is sensi- 
tive only to the macroscopic equation of state for the ISM. 
Whether or not the effective equation of state we have adopted 
(i.e. Figure^ for the ISM is appropriate in detail is uncertain, 
but this question is largely independent of the exact proper- 
ties of star-forming gas. In particular, we expect our results to 
generalize to other models for the ISM with similar effective 
equations of state. 

Having identified one exponential disk galaxy in our simu- 
lation, future work is needed to extend our analysis to other 
galaxies. For example, a good statistical sample of simulated 
disks is needed before we can characterize the distribution 
of disk sizes and their morphological types, as well as the 
photometric and kinematic properties of the galaxies. This 
work will require larger simulations of substantial dynamic 
range in order to obtain adequate numbers of simulated 
disk galaxies, but we find the initial successes described 
here and recently by other authors highly encouraging. A 
comparison with the wealth of observational data on the 
structural parameters of galaxies in the local universe will 
then show whether hydrodynamical simulations are finally 
producing realistic disk galaxies. 
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